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ABSTRACT 


We  are  developing  regional-phase  (Pn,  Pg,  Sn,  Lg)  tomographic  attenuation  models  for  Eurasia.  The  models  will  be 
integrated  into  the  National  Nuclear  Security  Administration  (NNSA)  Knowledge  Base  and  used  in  the  Magnitude 
and  Distance  Amplitude  Correction  (MDAC)  station  calibration  for  the  development  of  regional  seismic 
discriminants.  Our  current  focus  is  on  Pn,  an  extremely  important  phase  in  seismic  event  identification. 


Accurately  accounting  for  regional-phase  geometric  spreading  is  critical  for  the  development  of  useful  attenuation 
models.  It  is  particularly  important  for  Pn  and  Sn  waves  because  the  propagation  mode  of  these  waves  makes  them 
more  susceptible  to  upper  mantle  velocity  structures  and  the  Earth’s  sphericity,  which  in  turn  causes  the  geometric 
spreading  of  Pn  (and  Sn)  to  be  dependent  on  frequency  as  well  as  on  range  in  a  complicated  way.  We  conduct 
numerical  simulations  to  quantify  P„  and  S„  geometric  spreading  in  a  spherical  Earth  model  with  constant  mantle 
velocities.  Based  on  our  simulation  results,  we  present  new  P„  and  S„  geometric  spreading  models  in  the  forms 
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where  i  =  1,  2  or  3.  r  is  epicentral  distance;/is  frequency;  n,  =  1  km  and  f  ,  =  I  Hz.  We  derive  values  of  coefficients 
ntj  by  fitting  the  model  to  computed  P„  and  S„  amplitudes  for  a  spherical  Earth  model  having  a  40-km-thick  crust, 
generic  values  of  P  and  S  velocities  and  a  constant-velocity  uppermost  mantle. 

We  apply  the  new  spreading  model  to  observed  data  in  Eurasia  to  estimate  average  P„  attenuation,  obtaining  more 
reasonable  results  compared  to  using  the  standard  power-law  model.  Our  new  P„  and  S„  geometric-spreading  models 
provide  generally  applicable  reference  behavior  for  spherical  Earth  models  with  constant  uppermost-mantle 
velocities.  Since  our  region  of  interest  encompasses  most  of  Eurasia  from  the  Equator  to  the  North  Pole,  simple 
regular  gridding  methods  such  as  dividing  the  region  with  latitude  and  longitude  lines  would  result  in  cells  with 
drastically  different  cell  sizes  for  the  tomographic  inversion.  To  avoid  this  problem  and  to  optimize  the  resolution  of 
the  tomographic  model  based  on  data  distribution,  we  have  implemented  several  gridding  schemes.  They  include 
equal-cell-size  gridding  and  variable-cell-size  gridding  based  on  certain  criteria  such  as  the  number  of  path  hits. 
These  gridding  schemes  will  be  tested  in  future  tomographic  inversions  to  find  an  optimum  gridding  method. 

To  prepare  for  the  attenuation  tomographic  inversion,  we  have  been  collecting  and  measuring  regional-phase 
amplitude  data  from  the  Incorporated  Research  Institute  for  Seismology  Data  Management  Center  and  from  the  Los 
Alamos  National  Laboratory  Ground-Based  Nuclear  Explosion  Monitoring  Research  and  Engineering  program 
database.  Initial  Pn  amplitude  measurements  exhibit  similar  decay  behavior  as  what  we  see  in  the  synthetic  data.  To 
improve  our  data  coverage,  we  are  conducting  further  data  collection,  phase  picking  and  amplitude  measurement. 
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OBJECTIVE 


The  objective  of  this  project  is  to  develop  1-Hz,  two-dimensional,  regional-phase  {Pn,  Pg,  Sn  and  Lg)  tomographic 
attenuation  models  for  Eurasia.  The  models  will  be  used  in  MDAC  for  improved  event  identification. 

RESEARCH  ACCOMPLISHED 


Pn  geometric  spreading 

Accurately  accounting  for  geometric  spreading  is  critical  for  the  development  of  meaningful  regional-phase 
attenuation  models.  This  is  particularly  true  for  Pn  and  Sn  waves  because  the  nature  of  their  wave  propagation 
renders  them  acutely  sensitive  to  upper  mantle  velocity  structure  and  the  Earth’s  sphericity.  Even  simple  one¬ 
dimensional  (ID)  velocity  models  can  produce  geometric  spreading  of  Pn  and  Sn  that  is  strongly  dependent  on 
frequency  and  range  (e.g.,  Sereno  and  Given,  1990).  If  frequency  dependence  of  the  geometric  spreading  actually 
occurs  and  is  neglected,  the  attenuation  model  will  acquire  incorrect  frequency  dependence.  Similar  arrival  times  of 
Pg  and  Pn  phases  and  Pn  and  P  phases  at  their  respective  crossover  distances  result  in  rapidly  changing  P- wave 
amplitudes,  difficulty  in  phase  isolation  and  identification,  and  uncertainty  in  appropriate  specification  of  the 
propagation  path  and  geometric  spreading  at  these  distances.  Lateral  variation  of  Moho  topography  and  upper- 
mantle  lid  velocity  and  fine  scale  heterogeneity  of  the  lower  crust  and/or  mantle  lid  further  introduce  2D  and  3D 
complexities  into  Pn  and  Sn  spreading. 


We  conduct  most  of  our  simulations  using  the  reflectivity  method.  Results  from  reflectivity  calculations  are 
compared  with  results  calculated  with  a  2D  finite-difference  code  (Xie  and  Lay,  1994)  and  a  2.5D  axisymmetric 
spherical  finite-difference  code,  SHaxi,  to  confirm  that  the  earth  flattening  transform  (EFT)  and  layer  discretization 
required  by  the  reflectivity  method  do  not  produce  numerical  artifacts.  The  reflectivity  method  generates  complete 
synthetic  seismograms  within  a  specified  slowness  range  for  ID,  plane-layered  velocity  models.  In  order  to  use  the 
reflectivity  method  for  a  spherically  symmetric  Earth  model,  we  apply  the  EFT  to  transform  the  spherical  Earth 
model  to  a  plane  Earth  model.  Transformations  of  velocity  v  and  depth  z  are  (Chapman,  1973;  Muller,  1977) 
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where  R  is  the  radius  of  the  Earth.  Subscript  r  designates  values  in  the  spherical  (radially  symmetric)  model  and 
subscript/ designates  values  in  the  plane  (flat)  model.  The  density  p  transformation  is 
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which  is  not  unique  since  m  can  take  any  value  between  -5  and  1.  For  regional  body-waves,  the  choice  of  m  is  not 
critical  (Muller,  1977).  We  choose  m  =  -1  for  P/SV  simulations  (Muller,  1977)  and  m  =  -5  for  SH  simulations 
(Chapman,  1973).  We  experimented  with  different  values  of  m  and  the  results  were  basically  unchanged.  Finally  the 
transformation  of  amplitudes  calculated  from  plane-model  simulations  back  to  corresponding  amplitudes  in  the 
spherical  model  is 
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where  Zr  is  the  depth  of  the  source  in  the  spherical  model  and  A  is  epicentral  distance  in  radians. 


We  approximate  the  velocity  gradient  resulting  from  the  EFT  (Equations  1  and  2)  with  homogeneous  layers  in  the 
plane  Earth  model,  as  is  required  by  the  reflectivity  method.  The  thickness  of  these  layers  affects  the  accuracy  of  the 
approximation,  with  thinner  layers  yielding  more  accurate  results.  We  set  the  thickness  of  these  layers  to  be  about 
0.4  of  the  minimum  wavelength  of  the  waves  to  be  modeled,  which  appears  to  be  more  than  adequate.  Further 
reducing  the  ratio  (e.g.,  from  0.4  to  0.2  of  the  minimum  wavelength)  does  not  alter  the  results  appreciably.  The  total 
thickness  of  the  gradient  zone  is  set  to  be  more  than  100  km  larger  than  the  maximum  penetration  depth  of  the  direct 
wave  in  a  homogeneous  spherical  model  recorded  at  the  longest  epicentral  distance  considered.  This  thickness 
guarantees  that  no  P„  or  S„  waves  observed  within  the  distance  range  of  interest  are  affected  by  the  lower  boundary 
of  the  gradient  zone.  Below  the  gradient  zone,  the  velocity  is  constant. 
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We  use  the  same  generic  spherical  Earth  model  considered  by  Sereno  and  Given  (1990)  as  the  base  model  for  our 
simulations  and  use  the  synthetics  from  the  simulation  to  derive  parameters  of  the  P„  and  S„  geometric-spreading 
models  that  we  develop.  The  reference  Earth  model  consists  of  a  40-km-thick  outer  layer,  representative  of  an 
average  continental  crust,  with  a  constant-velocity  mantle  underneath  (Figure  1).  The  model  has  no  anelastic 
attenuation.  The  simplicity  of  this  model  allows  us  to  isolate  the  effects  of  the  sphericity  on  P„  and  S„  geometric 
spreading.  We  use  an  isotropic  point  source  for  P„  simulations.  The  source  for  S„(SH)  simulations  is  a  fundamental 
fault  vertical  strike-slip  source  and  the  source  for  Sn(SV)  simulations  is  a  dip-slip  source.  For  all  source  types  in  our 
main  calculations,  a  delta  function  is  used  as  the  source  time  function;  source  depth  is  15  km  and  source  strength  is 
10 15  N  m.  Three-component  synthetic  ground  displacements  are  computed  at  33  locations  distributed  log-evenly 
along  a  linear  profile  from  200  km  (1.8°)  to  2500  km  (22.5°).  The  Nyquist  frequency  of  the  seismograms  is  20  Flz. 


Figure  1.  Base  Earth  model  used  for  P„  and  S„  simulations  and  the  development  of  new  P„  and  S„  geometric¬ 
spreading  models.  Quality  factor  Q  is  infinite  throughout  the  model. 

We  cut  P„  and  S„  portions  of  the  synthetic  seismograms  using  fixed-velocity  windows.  The  velocities  that  we  use  to 
define  the  widths  ofP„  windows  are  7.6  km/sec  and  8.2  km/sec  and  those  for  S„  windows  are  4.0  km/sec  and  4.7 
km/sec  (Hartse,  et  ah,  1997).  The  windows  are  centered  at  the  peaks  of  the  phases.  We  also  tested  a  fixed-window- 
width  method  and  the  results  remained  essentially  the  same.  We  window  P„  and  Sn(SV)  from  vertical-component 
seismograms  and  S„(SH)  from  transverse-component  seismograms.  After  P„  and  S„  seismograms  are  windowed,  we 
taper  the  seismograms  with  small  tapers  (between  2%  to  20%  depending  on  the  length  of  the  signal  relative  to  the 
window  length)  and  Fourier  transform  the  seismograms  to  obtain  the  amplitude  spectra.  We  make  spectral- 
amplitude  measurements  at  100  frequencies  log-evenly  distributed  between  0.75  and  13  Flz.  Amplitude  at  each 

frequency  f  is  calculated  by  taking  the  average  of  the  amplitudes  between  frequencies  fl  -\f2  and  a/2 f. 

To  accurately  assess  the  geometric  spreading  of  seismic  phases,  the  propagation  medium  used  for  the  simulation 
should  have  no  attenuation.  Flowever,  in  order  to  avoid  a  computational  singularity,  the  reflectivity  method  requires 
a  nonzero  amount  of  attenuation  for  the  medium.  We  take  an  asymptotic  approach  similar  to  that  used  by  Yang 
(2002)  to  derive  P„  and  S„  amplitudes  for  an  elastic  model  without  attenuation  from  amplitudes  calculated  for  a 
group  of  anelastic  models.  We  first  make  20  calculations  for  models  that  have  attenuation  quality  factor  Q  log- 
linearly  increasing  from  10,000  to  100,000.  For  each  calculation,  a  single  Q  is  used  for  both  P  and  S  waves  and  for 
all  parts  of  the  model.  Amplitudes  at  each  frequency  and  each  epicentral  distance  from  these  calculations  are  then  fit 
by  a  quadratic  polynomial  as  a  function  of  HQ.  The  limit  of  the  polynomial  as  Q  approaches  infinity  is  taken  as  the 
amplitude  at  that  frequency  and  distance  for  the  elastic  model. 

P„  Modeling  Results 

Figure  2  plots  the  vertical  synthetic  P„  seismograms  from  the  base-model  simulation  at  selected  epicentral  distances. 
Q  used  in  this  simulation  is  100,000.  The  seismograms  are  low-pass  filtered  below  10  Flz  to  suppress  numerical 
noise  near  the  Nyquist  frequency.  The  figure  reveals  several  interesting  characteristics  of  P„  traveling  in  a  spherical 
Earth  model  with  constant  mantle  velocities.  Due  to  the  sphericity,  the  apparent  P„  velocity  is  not  constant,  but 
varies  with  epicentral  distance.  As  is  predicted  by  theory  (e.g.,  Cerveny  and  Ravindra,  1971),  the  pulse  shape  of  P„ 
evolves  from  that  of  the  impulse  source  at  distances  close  to  the  critical  distance  (about  0.8°  for  the  base  model  and 
a  15-km-deep  source)  to  the  shape  of  a  far-field  body  wave,  which  is  the  time  derivative  of  the  source  pulse,  at 
farther  distances.  The  amplitude  of  the  phase  changes  in  a  complex  manner,  first  decreasing  and  then  increasing, 
within  this  distance  range.  At  about  10°  to  12°,  the  first  pulse  separates  from  the  rest  of  the  Pn  wave  packet  and 
somewhere  between  16°  and  19°,  a  second  pulse  separates. 
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Figure  2.  Synthetic  P„  seismograms  from  base-model  calculations.  The  seismograms  are  filtered  below  10  Hz. 

Travel  time  is  reduced  by  8.2  km/sec.  Only  every  other  trace  calculated  is  plotted  to  enhance 
clarity,  r  is  epicentral  distance  in  kilometers. 


For  high-frequency  P„  at  distances  away  from  the  critical  distance,  Cerveny  and  Ravindra  (1971)  offer  a  detailed 
decription  of  the  signal  behavior  from  ray  theory.  Although  Cerveny  and  Ravindra  (1971)  describe  the  phenomenon 
for  a  plane-layered  Earth  model  with  positive  and  constant  velocity  gradient  in  the  mantle,  their  description  and 
conclusions  are  applicable  to  the  spherical-Earth-model  situation  as  well  since  the  spherical  model  can  be  mapped, 
through  the  EFT,  into  a  plane-layered  model  with  an  approximately  constant  velocity  gradient  in  the  uppermost 
mantle.  Following  Cerveny  and  Ravindra  (1971),  the  P„  phase  at  distances  between  about  5°  and  10°  in  Figure  2  can 
be  thought  of  as  the  superposition  of  individual  waves  reflected  n  times  (/?  =  0,  1,2,  . . .)  from  the  underside  of  the 
Moho.  The  superposed  wave  is  termed  interference  head  waves  by  Cerveny  and  Ravindra  (1971)  and  is  likened  to 
the  “whispering  gallery”  phenomenon  by  Menke  and  Richards  (1980).  As  distance  increases,  individual  components 
of  the  interference  head  wave  start  to  separate  from  the  wave  packet  due  to  their  increasingly  shorter  path  lengths 
compared  with  path  lengths  of  the  remaining  waves  in  the  wave  packet.  The  first  wave  to  separate  is  the  wave  that 
has  no  reflection  at  the  Moho  (the  direct  or  diving  wave).  This  is  evidenced  as  the  separation  of  the  first  pulse  in 
Figure  2.  The  second  separated  pulse  in  the  figure  is  the  wave  that  is  reflected  once  from  the  Moho.  From  ray 
theory,  the  epicentral  distance  at  which  the  /v-timc  reflected  wave  separates  from  the  interference  head-wave  packet 
is  (Cerveny  and  Ravindra,  1971,  Eq.  6.4) 
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where  H  is  crustal  thickness;  d  is  source  depth;  vc  is  P-wave  velocity  of  the  crust;  vm  is  P- wave  velocity  at  the  top  of 
the  mantle;  T  is  pulse  width  of  the  wave  and  g  is  velocity  gradient  ( dv/dz )  at  the  top  of  the  mantle  resulting  from  the 
EFT.  From  Equation  4  and  for  the  base  model,  a  wave  with  a  pulse  width  of  about  0.1  second  (10-Hz)  that 
undergoes  no  reflection  at  the  Moho  will  separate  at  about  9.8°  and  the  wave  that  has  one  reflection  will  separate  at 
about  16.6°.  These  predictions  are  consistent  with  the  synthetics  in  Figure  2. 


Figure  3  shows  the  amplitude  spectra  of  synthetic  Pn  at  the  same  epicentral  distances  as  those  in  Figure  2  from  the  Q 
=  100,000  computation.  The  figure  illustrates  the  evolution  ofP„  spectrum  from  proportional  to  the  source  spectrum 
to  proportional  to  the  time  derivative  of  the  source  spectrum. 
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Figure  3.  Amplitude  spectra  of  P„  at  the  same  distances  as  those  in  Figure  2.  Some  of  the  distances  are 
marked  on  the  left  of  the  corresponding  spectra. 


As  was  described  above,  we  use  an  asymptotic  method  to  derive  spectral  amplitudes  for  an  elastic  model  from 
amplitudes  obtained  using  anelastic  simulations.  Figure  4  gives  an  illustration  of  the  method.  Plotted  in  the  figure 
are  P„  amplitudes  at  different  frequencies  measured  from  calculations  using  different  quality  factors.  The  amplitudes 
are  computed  at  22.5°  epicentral  distance  and  are  normalized  by  the  maximum  amplitude  in  the  figure.  Quadratic 
polynomial  fits  of  the  amplitudes  are  plotted  as  solid  lines.  The  polynomial  fits  are  almost  perfect,  indicating  that 
our  approach  is  appropriate.  Amplitudes  at  other  distances  are  fit  as  well  as  those  shown  in  Figure  4. 
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Figure  4.  An  illustration  of  how  elastic  results  are  obtained  from  anelastic  calculations.  Circles  are 

amplitudes  at  22.5°  epicentral  distance  for  different  frequencies  plotted  against  Q  used  in  their 
calculations.  The  amplitudes  are  normalized  by  the  maximum  amplitude  in  the  plot,  which  is  the 
10.03-Hz  amplitude  from  the  Q  =  100,000  calculation.  Solid  lines  are  corresponding  quadratic 
polynomial  fits  of  the  amplitudes.  Numbers  on  the  right  are  the  limits  of  normalized  amplitudes  at 
corresponding  frequencies  as  Q  approaches  infinity. 


To  visualize  the  P„  amplitude  decay  in  a  spherical  Earth  model,  we  plot  10-FIz  P„  amplitudes  for  the  base  model  in 
Figure  5.  We  extend  the  epicentral-distance  range  to  between  135  km  (1.2°)  and  8000  km  (71.9°)  for  this  particular 
simulation  in  order  to  better  depict  the  evolution  of  Pn  waves.  Amplitudes  at  distances  beyond  about  20°  are 
measured  from  the  direct  wave  that  has  been  completely  separated  from  the  interference  head  waves.  The 
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amplitudes  are  corrected  for  the  free-surface  effect,  which  is  only  important  at  teleseismic  distances.  Also  plotted  in 
the  figure  are  the  amplitude  decay  of  a  conical  head  wave  in  a  plane  one-layer-over-half-space  model  (Aki  and 
Richards,  2002;  Eq.  6.26)  and  the  amplitude  decay  of  infinite-frequency  direct  wave  in  a  spherical  Earth  model  from 
ray  tracing.  At  distances  close  to  the  critical  distance,  P„  geometric  spreading  behaves  like  that  of  a  conical  head 
wave.  As  distance  increases,  P„  spreading  starts  to  deviate  from  that  of  the  head  wave  and  at  about  5°,  P„  amplitudes 
begin  to  increase.  As  was  mentioned  before,  10-Hz  direct-wave  energy  would  separate  from  the  rest  of  the 
interference  head  wave  at  about  10°.  It  seems  from  the  figure  that  this  separation  is  manifested  in  a  change  in  the 
smoothness  of  the  P„  amplitude  variation  followed  by  a  reduced  rate  of  amplitude  increase.  In  the  range  beyond  the 
critical  distance  and  before  the  direct-wave  separation,  P„  evolves  from  a  wave  similar  to  a  conical  head  wave  to  the 
interference  head  wave,  which  is  a  superposition  of  multiple  waves  reflected  from  the  Moho.  As  the  epicentral 
distance  approaches  teleseismic  distances,  the  direct-wave  spreading  approaches  that  of  the  infinite-frequency  wave 
from  ray  tracing  results,  as  is  expected.  The  direct  wave  dominates  the  whole  P„  wave  packet  at  long  distances.  We 
do  not  see  significant  difference  between  spectral  amplitudes  obtained  by  windowing  the  whole  P„  wave  packet  and 
those  obtained  by  just  windowing  the  direct  wave  after  its  separation  from  the  packet.  This  is  consistent  with 
theoretical  predictions  (Cerveny  and  Ravindra,  1971). 


Figure  5. 10-Hz  synthetic  P„  amplitude  decay  in  a  spherical  Earth  model  with  constant  mantle  velocities.  The 
solid  line  depicts  the  theoretical  amplitude  decay  of  a  conical  head  wave  in  a  plane  one-layer-over¬ 
half-space  Earth  model.  The  dashed  line  is  the  amplitude  decay  of  infinite-frequency  direct  wave  in 
a  spherical  homogeneous  Earth  model  from  ray-tracing  calculations. 


P„  geometric  spreading  in  a  spherical  Earth  model  is  not  only  different  from  that  of  a  head  wave  as  is  shown  in 
Figure  5,  but  also  frequency  dependent.  Figure  6  shows  the  Pn  amplitude-variation  surface  as  a  function  of  distance 
and  frequency  for  the  base  model.  The  strong  frequency  dependence  of  the  amplitudes  is  apparent.  Amplitudes  at 
higher  frequencies  are  affected  more  by  the  sphericity  than  are  lower-frequency  amplitudes.  The  separation  distance 
of  the  direct  wave  from  the  interference  head  waves  becomes  shorter  as  frequency  becomes  higher  (Equation  4). 
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Figure  6.  Synthetic  P„  amplitudes  as  a  function  of  epicentral  distance  and  frequency. 


A  New  P„  Geometric-Spreading  Representation  for  a  Spherical  Earth  Model 

Figures  5  and  6  illustrate  that  a  frequency-independent,  power-law  model  cannot  accurately  represent  P„  geometric 
spreading  in  a  spherical  Earth  model.  Such  a  representation  would  plot  as  a  straight  line  in  Figure  5,  which  is  clearly 
inappropriate  for  modeling  P„  geometric  spreading  over  a  wide  distance  range.  In  addition,  a  power-law  model  with 
constant  exponent  does  not  take  into  account  the  frequency  dependence  of  P„  spreading  shown  in  Figure  6.  Based  on 
the  P„  amplitude-decay  behavior  shown  in  Figures  5  and  6,  we  propose  a  new  empirical  P„  geometric-spreading 
model  that  fits  the  synthetic  data  much  better  and  that  also  results  in  more  reasonable  anelastic-attenuation  estimates 
from  observed  data,  as  we  will  discuss  in  more  detail  in  the  next  section. 


(5) 

(6) 

(7) 


In  Equation  5,  A  is  a  frequency-dependent  scaling  factor;  M0  is  source  moment;  R  is  source  radiation  pattern;  Q  is  Pn 
quality  factor;  v  is  P„  velocity;  S  is  receiver  site  response;  r  is  epicentral  distance;  <9  is  azimuth  angle  and/is 
frequency.  r0  and  /«  are  included  in  Equations  6  and  7  in  order  for  the  new  model  to  have  the  same  dimension  as 
standard  power-law  models  (e.g.,  Street  et  al.,  1975;  Sereno  et  al.,  1988).  The  main  differences  between  the  new 
geometric-spreading  model  (Equations  6  and  7)  and  the  standard  frequency-independent  power-law  model  are  the 
addition  of  the  first  term  in  the  exponent  and  the  frequency  dependence  of  parameters  In  the  logarithm  domain, 
the  new  model  is  a  quadratic  function  of  log-distance,  whereas  the  power-law  model  is  linear.  The  reason  for 
choosing  a  log-quadratic  function  is  to  keep  the  parameterization  as  simple  as  possible  while  providing  a  good  fit  to 
the  synthetics.  The  adoption  of  a  quadratic  functional  form  for  «,■  (Equation  7)  is  based  on  the  behavior  of  n,  versus 
frequency  obtained  by  fitting  Equation  6  to  synthetic  Pn  amplitudes  at  individual  frequencies. 
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If  we  take  common  logarithm  of  Equation  6,  substitute  Equation  7  into  the  result  and  let  r0  and  fa  equal  one,  we 
obtain 

log  [G(r,f)\  =  «„(log/)2(logr)2  +  nl2  log/(logr)2  +  «l3(logr)2 

-«2i  (log/)2logr-«22  log/log  r-n23  log  r  ,  (8) 

+»3i(log/)2  +  n32  log/  +  «33 

where  r  is  in  kilometers  and/is  in  hertz.  To  derive  coefficients  tty,  we  fit  Equation  8  to  synthetic  P„  amplitudes 
shown  in  Figure  6  in  a  least-squares  sense.  P„  amplitudes  are  corrected  for  M0  used  in  the  simulation  (1015  N  m)  and 
K  before  the  fitting.  Since  the  source  that  we  use  in  the  simulation  has  a  flat  spectrum,  K  is  frequency  independent. 
We  use  K  =  (4  7tpv)~l  (Denny  and  Johnson,  1991),  where  p  is  density  and  v  is  P-wave  velocity  of  the  source  region. 
Source  radiation  and  site  response  are  unity.  We  use  P„  amplitudes  at  epicentral  distances  beyond  300  km  (2.7°)  and 
before  the  start  of  the  direct- wave  separation  to  fit  the  model.  We  use  300  km  as  the  lower  distance  limit  because 
reliable  P„  observations  are  typically  made  at  some  distances  beyond  the  Pg  crossover  distance  (-200  km).  The 
choice  of  300  km  is  also  to  avoid  possible  long-period  numeric-noise  contamination  at  short  distances,  as  is 
indicated  in  Figure  6.  The  upper  distance  limits  are  based  on  the  observation  that  within  these  distances,  P„  is  the 
result  of  the  interference  of  all  of  its  components  including  the  direct  wave.  At  larger  distances,  the  direct  wave 
separates  from  the  rest  of  the  wave  packet  and  the  characteristics  of  P„  become  different.  The  P„  amplitude  decay 
within  the  defined  distance  range  also  has  a  smooth  pattern  and  thus  is  easier  to  fit  by  a  simple  mathematical  model. 
The  upper  distance  limits  vary  from  7.3°  to  17.3°  for  the  frequency  range  between  13  and  0.75  Hz.  Within  the 
specified  distance  limits,  the  new  spreading  model  is  applicable.  Since  at  about  15°,  P„  in  the  real  Earth  is  overtaken 
by  upper-mantle  triplications  resulting  from  reflections  and  refractions  at  4 1 0-km  and  660-km  discontinuities  and  is 
no  longer  the  first  arrival,  P„  is  usually  used  within  the  distance  range  where  the  new  spreading  model  is  valid  for 
frequencies  below  about  2  Hz.  For  higher  frequencies,  the  range-of-applicability  of  the  new  model  is  shorter,  but 
observationally  high-frequency  signals  are  generally  only  detectable  above  the  noise  level  at  shorter  distances. 
Coefficients  rty  (i  =  1,  2,  3;j  =  1,  2,  3)  from  the  fitting  are  listed  in  Table  1.  The  inclusion  of  r0  and  in  the  model 
also  guarantees  that  even  though  the  values  of  the  coefficients  are  derived  using  Equation  8  with  r  in  kilometers  and 
/in  hertz,  they  are  valid  for  r  and/in  any  units  as  long  as  r0  and/)  are  converted  accordingly. 


Table  1  Coefficients  of  the  new  P„  geometric-spreading  model 


«11 

«12 

«13 

«21 

«22 

«23 

«31 

«32 

/J33 

-0.217 

1.79 

3.16 

-1.94 

8.43 

18.6 

-3.39 

9.94 

20.7 

Application  to  Observed  Data 

The  key  value  of  any  mathematical  model  of  the  physical  world  is  for  the  model  to  be  able  to  provide  physically 
reasonable  descriptions  of  observed  data.  To  test  the  validity  and  usefulness  of  the  new  Pn  geometric-spreading 
model  and  the  associated  choice  of  a  constant-velocity  mantle  lid  structure,  we  correct  a  set  of  observed  P„  spectral 
amplitudes  for  geometric  spreading  with  the  new  model  and  estimate  the  average  medium  attenuation.  We  then 
compare  the  results  with  those  published  in  the  literature. 

We  represent  observed  P„  amplitudes  by  Equation  5.  For  the  purpose  of  testing  the  new  Pn  geometric-spreading 
model,  we  simplify  Equation  5  by  assuming  that  site  response  is  unity  for  all  stations  and  source  radiation  patterns 
can  be  ignored.  We  presume  that  errors  introduced  by  these  simplifications  are  random  and  should  not  affect 
average-attenuation  estimates  systematically.  With  known  or  estimated  source  moments,  an  assumed  scaling  factor 
K  and  a  P„  geometric-spreading  model,  we  can  estimate  the  average  attenuation  quality  factor  at  each  frequency  by 
least-squares  fitting  the  logarithm  of  source  and  geometric-spreading  corrected  spectral  amplitudes  as  a  function  of 
epicentral  distance. 

The  observed  P„  amplitudes  are  measured  on  vertical-component  ground-displacement  data  recorded  by  stations  in 
and  around  China  and  in  southern  Europe  for  events  in  the  same  region.  The  same  windowing  method  as  the  method 
we  employ  to  measure  the  synthetic  P„  amplitudes  is  used.  Analyst  picks  reported  in  global  catalogs  (ISC,  EDR, 
REB,  EHB,  etc.)  are  used  to  center  the  P„  windows.  We  derive  source  moments  from  body-wave  magnitudes  (////,). 
We  only  use  amplitudes  from  events  with  nib  equal  to  or  smaller  than  6  to  avoid  magnitude  saturation.  We  use  a 
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simplified  version  of  the  scaling  factor  K  expressed  as  K(f)  =  {4^ov3[l+(///^)2]}_1  using  crustal  P-  wave  velocity  and 
density  of  the  base  Earth  model  as  v  and  p.  Source  corner  frequency  /,  is  calculated  from  source  moment  using  the 
relationship  logM0  =  1 7.08-3.24log/,.  derived  by  Xie  and  Patton  (1999)  from  P„  amplitude  data  recorded  in  central 
Asia.  For  comparison,  we  use  both  the  new  Pn  geometric-spreading  model  and  the  power-law  model  with  two 
different  exponents,  -1.1  (Walter  and  Taylor,  2002)  and  -1.3,  in  the  attenuation  estimation.  When  the  new  spreading 
model  is  used,  we  limit  the  epicentral  distances  of  the  amplitudes  used  in  the  estimation  to  within  the  distance  range 
where  the  model  is  valid.  For  power-law  model  corrections,  we  use  amplitudes  between  300  km  (2.7°)  and  1668  km 
(15°).  Attenuation  is  estimated  at  0.75,  1.0,  2.0,  4.0  and  6.0  Flz.  The  new  spreading  model  is  used  in  the  correction. 
Although  the  amplitudes  show  a  large  scatter,  a  linear  decay  trend  due  to  realistic  values  of  medium  attenuation  is 
discernable. 

Table  2  lists  estimated  average  quality  factors  using  different  geometric-spreading  corrections  from  observed 
P„  amplitudes.  Using  the  power-law  spreading  model  results  in  some  negative  values  at  low  frequencies.  At  high 
frequencies,  the  power-law  model  yields  estimates  that  range  from  over  1000  to  over  5000.  On  the  other  hand, 

Q  estimates  using  the  new  spreading  model  are  positive  and  below  700  at  all  frequencies.  These  values  can  be 
compared  with  published  P- wave  quality  factors  in  the  mantle  lid,  as  we  discuss  below. 


Table  2  Average  P„  quality  factor  estimated  using  different  geometric-spreading  models 


data  frequency  (Flz) 

0.75 

1.0 

2.0 

4.0 

6.0 

new  model 

440 

338 

312 

557 

678 

power-law  model  1* 

-1669 

9241 

1206 

2204 

3450 

power-law  model  2** 

-734 

-2153 

1850 

3260 

5246 

model  with  exponent  of  -1.1 


model  with  exponent  of  -1.3 


Using  theory,  observed  body-wave  spectra  and  waveform  modeling,  Lundquist  and  Cormier  (1980)  derive  generic 
absorption-band  P- wave  Q  models  for  the  mantle.  The  Q  values  of  these  models  range  from  about  100  to  500  for 
frequencies  between  0.7  Flz  and  6.0  Flz  in  the  depth  range  of  45  to  200  km.  In  their  paper,  Lundquist  and  Cormier 
(1980)  also  cites  results  of  some  other  Q  studies  that  use  free-oscillation  data,  long-period  surface  waves  and  high- 
frequency  (1-5  Flz)  body  waves.  The  frequency-independent  P- wave  Q  models  from  these  studies  have  values  from 
about  100  to  250  for  depths  between  50  and  150  km.  Der  et  al.  (1986)  construct  a  P- wave  Q  model  for  the  Eurasian 
Shield  using  a  large  set  of  teleseismic  body  waves.  Their  model  has  values  between  about  350  and  900  for 
frequencies  between  0.3  and  10  Flz  at  depths  between  100  and  200  km.  Above  100  km,  Q  values  increase  to  between 
600  and  1500  for  the  same  frequency  range. 

More  recently,  some  studies  make  direct  P„  Q  estimations.  Sereno  et  al.  (1988)  and  Sereno  (1990)  obtain  P„  Q 
models  for  Scandinavia  and  eastern  Kazakhstan  respectively  by  inverting  broadband  P„  spectra. 

The  0.75-  to  6-Flz  Pn  Q  values  that  they  estimated  are  between  283  and  768  for  Scandinavia  and  between  260  and 
735  for  eastern  Kazakhstan.  Although  Sereno  et  al.  (1988)  and  Sereno  (1990)  assume  a  power-law  P„  geometric¬ 
spreading  model  with  an  exponent  of -1.3,  their  P„  Q  estimates  are  more  in  line  with  the  average  P„  Q  estimates  that 
we  obtain  using  the  new  P„  spreading  model  than  with  those  from  power-law  model  corrections  (Table  2).  A 
possible  explanation  for  this  observation  is  that  the  majority  of  their  data  are  recorded  within  1000-km  epicentral 
distance.  At  short  distances,  the  power-law  spreading  model  has  a  gentler  slope  than  the  new  spreading  model  does 
and  therefore  would  yield  smaller  Q  estimates  from  short-distance  data.  Flowever,  for  a  broader  distance  range  such 
as  the  distance  range  that  our  dataset  covers,  the  power-law  model  yields  larger,  sometimes  negative,  Q  estimates 
because  of  the  steeper  slope  of  the  model  at  long  distances.  The  implication  is  that  if  a  power-law  P„  spreading 
model  with  a  specific  exponent  is  used,  it  will  only  be  applicable  as  an  approximation  in  a  limited  distance  range  and 
models  with  different  exponents  are  needed  for  different  distance  ranges.  Our  parameterization  remedies  this  failing. 
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Comparing  Q  values  in  Table  2  with  those  published  in  the  literature,  we  conclude  that  the  new  P„  geometric¬ 
spreading  model  yields  Q  estimates  that  are  generally  consistent  with  published  results  over  the  broad  region  of 
Eurasia.  The  Q  estimates  with  power-law  model  corrections,  on  the  other  hand,  have  values  that  are  either  negative 
or  seem  to  be  too  large.  It  should  be  noted  that  Q  values  estimated  in  this  fashion  represent  only  the  average  P„ 
attenuation  behavior  over  the  broad  region  of  Eurasia.  They  may  not  be  accurate  for  any  particular  region  due  to  the 
deviation  of  the  uppermost-mantle  velocity  structure,  including  any  non-zero  mantle  lid  velocity  gradient  or  lateral 
heterogeneity,  from  the  constant-velocity  structure  that  we  use  to  develop  the  geometric-spreading  model.  The  fact 
that  we  obtain  reasonable  Q  values  for  data  spanning  a  large  distance  range  demonstrates  that  our  reference  model  at 
least  provides  reasonable  average  behavior  without  any  peculiarities  such  as  negative  Q  values. 

S„  Simulations 

In  addition  to  simulating  P„  propagation  in  a  spherical  Earth  model,  we  also  simulate  S„  propagation  in  the  same 
base  model.  Except  for  different  source  types  and  different  slowness  integration  limits,  other  modeling  parameters  in 
the  Sn  simulation  are  kept  the  same  as  those  used  in  the  P„  simulation.  The  behavior  of  S„  in  a  spherical  Earth  model 
is  very  similar  to  the  behavior  ofP„  waves.  The  only  difference  is  that  the  separation  of  individual  waves  from  the 
interference  wave  packet  occurs  at  shorter  distances  for  S„.  This  difference  can  be  predicted  using  Equation  4, 
although  the  equation  was  originally  derived  only  for  P„  waves.  S„(SH)  spectral  amplitudes  also  form  a  surface  with 
a  shape  similar  to  that  of  the  P„  amplitude  surface  shown  in  Figure  6. 

Because  of  the  similarities  between  P„  and  S„  propagation  in  a  spherical  Earth  model  with  constant  mantle 
velocities,  we  propose  an  S„  geometric-spreading  model  that  has  the  same  functional  form  as  that  of  the  new  P„ 
spreading  model  (Equations  6  and  7).  We  derive  the  values  of  the  coefficients  tiy  by  fitting  synthetic  S„(SH) 
amplitudes.  S„(SV)  amplitudes  are  severely  contaminated  by  P- wave  energy  up  to  1000  km,  and  thus  are  not  suitable 
for  fitting.  Beyond  1000  km,  S„(SV)  and  S„(SH)  amplitudes  decay  similarly.  It  suggests  that  the  S„  spreading  model 
developed  by  fitting  S„(SH)  amplitudes  is  also  suitable  for  describing  S„(SV)  geometric  spreading. 

We  correct  synthetic  S„(SH)  amplitudes  for  source  moment,  source  radiation  pattern,  and  the  scaling  factor  K  before 
the  fitting.  We  use  K=  (\Kpv)\  but  now  v  is  5-wave  velocity  of  the  crust.  We  again  set  the  lower-distance  limit  to 
300  km  (2.7°)  for  amplitudes  used  in  the  fitting.  The  upper  distance  limits  are  from  6.2°  to  14.4°  for  frequencies 
from  13  Hz  to  0.75  Hz.  The  limits  set  the  distance  range  within  which  the  S„  spreading  model  is  valid.  Table  3  lists 
the  coefficients  n  j  for  the  S„  geometric-spreading  model  from  the  fitting. 


Table  3  Coefficients  of  the  new  S„(SH)  geometric-spreading  model 


«11 

n  12 

«13 

nn 

nn 

nn 

«31 

nn 

»33 

-0.347 

2.16 

3.54 

-2.69 

10.1 

20.4 

-4.38 

11.7 

23.1 

CONCLUSIONS  AND  RECOMMENDATIONS 


The  new  P„  and  S„  geometric-spreading  models  are  useful  in  common  situations  where  only  simple  velocity  models 
with  uppermost-mantle  structure  represented  as  constant-velocity  half  space  or  a  stack  of  constant  velocity  layers  are 
available.  If  the  mantle-lid  velocity  gradient  is  well  resolved  in  a  given  region,  which  is  sometimes  the  case, 
simulations  for  that  gradient  can  be  performed  to  obtain  appropriate  geometric-spreading  corrections.  Since  the 
effect  of  sphericity  is  equivalent  to  the  effect  of  a  positive  velocity  gradient  in  a  plane-layered  model,  we  anticipate 
that  the  functional  form  of  our  geometric-spreading  models  will  remain  the  same  for  Earth  models  in  which  an 
effective  (physical  plus  effects  of  the  sphericity)  positive  velocity  gradient  exists.  Only  the  coefficients  will  differ. 
Of  course,  any  Q  determinations  will  trade-off  directly  with  errors  due  to  having  the  wrong  velocity  model,  and  if 
one  can  determine  specific  structural  parameters  they  should  be  used  to  make  specific  geometric  spreading  terms. 

Quantifying  the  effects  of  specific  mantle-lid  velocity  gradients,  Moho  irregularity  and  lateral  velocity  heterogeneity 
in  the  uppermost  mantle  on  P„  and  S„  geometric  spreading  through  2D  and  3D  numerical  modeling  will  be  the 
subject  of  a  separate  study.  It  is  reasonable  to  assume  that  except  for  the  case  of  different  lid  velocity  gradients, 
many  effects  will  contribute  primarily  to  scatter  around  the  fundamental  behavior  of  our  geometric-spreading 
models. 
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